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SUMMARY 

A viscous parabolic subset of the incompressible Navier- 
Stokes equations, the quasi-cylindrical vortex equations, can be 
used to compute vortex flows as long as their stream surface 
angle remains small. At low values of a swirl parameter, 
computation presents no difficulties. At high swirl values, 
singularities are encountered which indicate a failure of the 
quasi-cylindrical approximation and are thought to be associated 
with physical axisymmetric vortex breakdown. 'Vortex breakdown' 
is normally followed by a 'vortex bubble', and a 'vortex jump' 
may occur farther downstream. 

It is found that the quasi-cylindrical vortex equations have 
an infinite number of discrete swirl-dependent singularities 
whose corresponding critical swirl values depend on velocity and 
circulation profile shapes. The behavior of flow (deceleration 
vs. acceleration on the vortex axis) is opposite on both sides 
of a singularity (i.e. at higher or lower swirl). At high swirls, 
the approach to singularities can only be avoided by application 
of specific external axial velocity and/or circulation gradients. 
The singularities appear to correspond to the critical swirl 
values of the equation of inviscid rotating flow. 

The first singularity is of particular significance as 
it seems to correspond to the often observed axisymmetric 
explosive breakdown. The physical significance of the higher 
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singularities is not quite clear at present. A second important 
value of the swirl parameter is S Q which separates flow which 
decays smoothly from vortex flow which eventually breaks down. 

Typical vortex flows (initially uniform axial flow, leading 
edge vortex, trailing vortex) are computed for a wide variety of 
initial swirl parameters using a method of weighted residuals 
which expresses the velocity and circulation approximations in 
terms of exponentials. The effects of external axial velocity 
and circulation gradients are investigated. 



I 


TABLE OF CONTENTS 

Page 

I. INTRODUCTION 1 

1.1 Physics of the Problem 1 

1.2 Swirl Parameter 8 

1.3 Numerical Approach 8 

1 . 4 Overview . H 

II. COMPUTATIONAL METHOD 13 

2.1 Governing Equations I 3 

2.2 Integral Relations 

2 . 3 Approximating and Weighting Functions .... 1® 

2.4 System of Ordinary Differential Equations . . 17 

2.5 Initial Profiles 18 

2.6 Accuracy and Convergence I 8 

III. SINGULARITIES AND CRITICAL SWIRL PARAMETERS ... 21 

3.1 Singularities of the Governing System .... 2 1 

3.2 Swirl Parameter 23 

3.3 Singularities for Initially Uniform 

Axial Flow 23 

3.4 Behavior of Axial Derivatives 23 

3.5 Effect of Velocity Profile 23 

3.6 Assessment of Breakdown Behavior 27 

IV. VORTEX COMPUTATIONS FOR CONSTANT EXTERNAL AXIAL 

VELOCITY AND CIRCULATION 29 

4.1 Initial Profiles and Distance to Failure . . 29 

Uniform initial axial flow 28 

Leading edge vortex ....... 31 

Trailing vortex 3 1 

Possibility of breakdown-^stable 

solutions of type 2, 3,4, 33 


v 


4.2 Velocity on the Axis 34 

Behavior of initially uniform axial flows 
of type 1-4 34 

Initially uniform axial flow of type 1 . . 36 

Leading edge vortex of type 1 . . 37 

Trailing vortex of type 1 38 

4.3 Velocity Profiles 38 

V. EFFECTS OF EXTERNAL VELOCITY AND CIRCULATION 

DISTRIBUTION 44 

5.1 Type 1 Vortex Flow 44 

5.2. Type 2 Vortex Flow 48 

5 . 3 Type 3 Vortex Flow 49 

5.4 Type 4 Vortex Flow 49 

5.5 Avoidance of Singularities 51 

VI. CONCLUSIONS 55 

APPENDIX 59 

REFERENCES ’ 6 3 


vi 



LIST OF SYMBOLS 


n 


n 


B 


C, 


D, 


E, 


B 


n ,k 
n ,k 


parameter in the velocity approximation 
parameter in the circulation approximation 


n ,k 
n ,k 


coefficients in the ordinary differential equa- 
ls tions (defined in the Appendix) 


= 2.71828183 


f k (V) 

g k Or> 

H 

j ln 

k 

K 

K 

e 

K e . 

e i 

N 


weighting functions 
weighting functions 
= VR (nondimensional parameter) 
zeros of the Bessel function 
= wr circulation 

= WR nondimensional circulation 
nondimensional external circulation 
initial K e 

number of parameters in the velocity and 
circulation approximations 


[-] 

[-] 


[-] 


[-] 

[-] 

(-] 

[-] 

[-] 

[l 2 /t] 

[-] 

[-1 

[-3 

(-3 


vii 


p pressure [F/L 2 ] 

p e external static pressure [F/L 2 ] 

p Q free-stream total pressure [F/L 2 ] 

P = p/pik, 2 /2) nondimensional pressure [-] 

P e nondimensional external static pressure [-] 

P Q nondimensional free-stream total pressure [-] 

nondimensional free-stream static pressure [-] 

r radial coordinate [L] 

r vortex core radius [L] 

c 

R = y/iie r/r c nondimensional radial coordinate [-] 

R c nondimensional vortex core radius [-] 

Re = lioTg/v - core Reynolds number [-] 

S swirl parameter S = Wq/u^ = W ( ,/U ax 

or s = (aw/aiO^/u^ [-] 

S(X) x-dependent swirl parameter as function of [-] 

velocity profile development 

S Cr values of S at which singularities exist [-] 

initial swirl parameter 

S D critical swirl parameter value (velocity profile- [-] 

dependent) separating vortex flow which decays 
smoothly (S<S 0 ) from vortex flow which eventually 
"breaks down" (S Q <S<Si) 

. . • singular swirl values (velocity profile- [-] 

dependent) where singularities of quasi- 
cylindrical vortex equations occur 

u axial velocity (x-direction) [L/T] 

u e external axial velocities [L/T] 

free-stream axial velocity [L/T] 

U = u/Ho nondimensional axial velocity [-] 

U g nondimensional external axial velocity [-] 

U ei initial Ug [-] 

v radial velocity (r-direction) [L/T] 

viii 



V 

w 


W C 

W 

W c 

x 

X 

x f 

X o 

Y 


a 

Y 

n 

v 

P 

a k 

4>(Y) 


» 


external radial velocity 

= yj Re (v/i^ ) nondimensional radial velocity 
swirl velocity (circumferential direction) 
swirl velocity at core radius r c 
= w/^ nondimensional swirl velocity 
= 

axial coordinate 

= x/r c nondimensional axial coordinate 

nondimensional position of computational failure 

x value at which velocity gradient is applied 

= R 2 /2 nondimensional radial coordinate 

exponent in velocity and circulation approximations 

exponent in the circulation profile 
-aY 

= e 

kinematic viscosity 
density 

weighting function exponent 

complete set of linearly independent functions 

d/dx 

d/dY 


Subscripts : 

k integer, 1 - k - N + 1 

n integer, 1 - n - N 

ax on the axis 


[L/T] 

[-] 

[L/T] 

[L/T] 

[-] 

[-] 

[L] 

[-] 

[-] 

[-] 

[-] 

[-] 

[-] 

[L 2 /T] 

[M/L 3 ] 

[-] 

[-] 

[-] 


ix 



I . INTRODUCTION 

1.1 Physics of the Problem 

The present report presents and discusses numerical 
solutions to a specific parabolic subset - applicable to quasi- 
cylindrical vortex flows - of the axisymmetric incompressible 
Navier-Stokes equations . This subset can be used as a reasonable 
approximation for many vortex flows occurring in technological 
applications and in nature (wing vortices, rotating tubes, 
swirling pipe flow, dust devils, tornados). Solutions of the 
quasi-cylindrical vortex equations are therefore of some interest 
in their own right, but of even greater interest are perhaps the 
conditions under which solutions cannot be obtained. In such 
cases the physics of the problem are different than assumed in 
the quasi-cylindrical approximation, and this subset becomes 
inadequate. When this happens, the physical counterpart of the 
flow under consideration most likely experiences a rapid expan- 
sion or contraction of the core to violate the condition of 
quasi-cylindricity . 

Observations of vortex flows at high swirl do indeed show 
rapid expansions and/or contractions under certain conditions. 

The observed closely-related phenomena have been lumped under 
the headings of 'vortex breakdown' and 'vortex bursting' , and 
several explanations have been advanced. Until recently a 
peculiarity of vortex flows - their extreme sensitivity to probe 
insertion - has prevented the gathering of reliable experimental 
data to verify these explanations. Only now, with the advent of 
probeless velocity measurement using the laser Doppler anemometer, 
can reliable data be collected, and theories be checked and 
assigned their proper place in the overall flow picture. It 
turns out - once more - that the conflicting theories are not 
really in conflict after all, but that they describe different 
phenomena which are indeed observed to often happen more or less 
under the same circumstances. 
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In order to put the results of this study into the proper 
physical perspective, an interpretation of 'vortex breakdown' 
observations will first be made partly based on recent laser 
anemometer studies which have lead to some unexpected observa- 
tions. The conclusions contrast some interpretations still to be 
found in the current literature. 

Earlier experimental and theoretical work in the area of 
concentrated vortex flows has been reviewed in Hall (1966a) , and 
a very comprehensive bibliography to 1967 is given in Timm 
(1967) . A review of vortex breakdown has recently been made by 
Hall (1972) , and a comprehensive review of confined vortex flows 
is due to Lewellen (1971) . Observations of wing vortex flows 
and vortex breakdown, especially on delta wings, are numerous 
!- references can be found in the reviews cited. With very few 
exceptions (especially Hummel 1965, McCormick et al. 1968) 
these studies are qualitative. It is expected that reliable 
quantitative data for wing vortex flows will become available 
in the near -future with the advent of two- and three-dimensional 
and rapid scanning laser Doppler anemometers . 

Much useful information on the phenomenon of vortex 
breakdown has come from the study of swirling flow in pipes (in 
particular Harvey 1962, Sarpkaya 1971a, 1971b, Orloff 1971, 
Orloff and Bossel 1971) . The dye studies of Harvey and 
Sarpkaya are qualitative only, and the conclusions of these 
authors have been somewhat revised by the quantitative optical 
velocity measurements of Orloff and Bossel. More information on 
the vortex breakdown phenomenon has come from the study of flow 
in a stationary cylinder with a rotating lid (Maxworthy 1967, 
Vogel 1968) . In the careful and comprehensive work of 
Vogel quantitative data on vortex breakdown are obtained from 
tracer studies. Swirling flows undergoing a rapid expansion 
have been studied by Gore and Ranz (1964) , Nissan and Bresan 
(1961), Potter et al . (1958), So (1967), Vonnegut (1954) and 

others . 


2 



Several conflicting explanations have been advanced to 
explain the sudden expansion of vortex cores , most commonly 
referred to as 'vortex breakdown' s the spiral instability theory 
of Ludwieg (1962, 1965), the hydraulic jump analogy of Benjamin 
(1962) , or, in a milder form, the concept of a stationary wave 
(Benjamin 1967, Leibovich 1968), and finally the explanation on 
the basis of stationary continuous solutions of the Navier-Stokes 
equations for swirling flow, or a proper subset thereof (Vaisey 
1956, Lavan and Fejer 1966, Bossel 1967, 1969, Orloff 1971, 

Orloff and Bossel 1971, Torrance and Kopecky 1971) . On the 
Other hand, the failure of the quasi-cylindrical subset of the 
Navier-Stokes equations (i.e. the failure of the quasi-cylindri- 
cal approximation in regions of rapidly expanding or contracting 
vortex flow) has been used by Gartshore (1963) , Hall (1966b) , 
Bossel (1967, 1971), Mager (1971) to predict the position of 
likely breakdown and the behavior of flow preceding it. 

Most experimental observations seem to support the view 
that the sudden expansion of vortex cores is most often an 
axisymmetric phenomenon and not the result of spiral instability. 
However, according to unpublished experiments (Ludwieg 1971, 
private communication) , spiral instability may indeed play a 
major role in the breakdown of delta wing vortices. Other 
instabilities have been identified by Sarpkaya (19 71a, b) , but the 
dominant role in what we know as the 'breakdown' -process appears 
to be a sudden axisymmetric expansion of the core. It is for 
this reason that the axisymmetric equations are studied in the 
present work. Before entering into the analysis, it will be 
useful to piece together a comprehensive picture of the axisym- 
metric 'vortex breakdown ' -process based on the most recent 
quantitative data obtained through optical velocity measurements 
by Orloff (1971) and Orloff and Bossel (1971) . It is found that 
the explanations and theories proposed so far all have their 
proper place and are not really competitive. The proper range 
of application will be identified. 
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Consider first a vortex flow with more or less cylindrical 
external stream surfaces and a high 'swirl' (some ratio of a 
representative circumferential (swirl) velocity to a representa- 
tive axial velocity) . An adverse pressure gradient is 
superimposed on the axis by the presence of an axisymmetric 
obstacle (Orloff and Bossel 1971) . The flow observed is then as 

shown in Fig. 1.1. The flow initially approaches the obstacle in a 
quasi-cylindrical manner. The velocity on the axis decreases 
until a stagnation point is reached on the axis and flow has to 
move outward from the axis near the stagnation point. In keeping 
with traditional notation, this particular process will be called 
'vortex breakdown' in the following. 

There is reversed (upstream) axial flow on the downstream 
side of the stagnation point. This reversed flow is part of 
the flow in the ('forced') 'vortex bubble' between stagnation 
point and obstacle, a closed region with practically no inter- 
change of fluid with the surrounding flow. The external flow 
passes over the bubble and over the rear of the obstacle. As 
it passes over the rear it undergoes a 'vortex jump' in the 
sense of the hydraulic jump analogy of Benjamin (1962) , with a 
substantial loss of axial momentum and an attendant expansion 
of stream surfaces near the axis . 

In Fig. 1.2 the stagnation point, and the rapid expansion . 
near the axis are caused by the effect of an external pressure 
gradient such as produced by a flow divergence. With a proper 
favorable pressure gradient the bubble will contract downstream 
and may even be more or less closed. Flows having such 
(seemingly) closed 'free' vortex bubbles can either be generated 
by proper shaping of an outer surface as in Fig. 1.2 (Bossel 
and Orloff 1971, unpublished) or by a self-induced process 
between vortex bubble and external flowfield as in the vortex 
tube flows of Harvey (1962) and Sarpkaya (1971a, b) . Sarpkaya 
has even succeeded in generating several vortex bubbles in a 
row. Just downstream of a vortex bubble a vortex jump may 
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occur, but does not appear to occur in many cases where the 
bubble is smooth or repeated. Following Benjamin (1967) and 
Leibovich (1968) these vortex bubbles can be viewed as solitary 
waves. However, their description appears to be more straight- 
forward and accessible to numerical quantitative analysis if 
direct use is made of the Navier-Stokes equations or a proper 
subset of these equations . 

Fig. 1.3 shows another . possibility which appears to be 
typical of delta wing vortex flows and some tornados. In this 
case the vortex bubble following the vortex breakdown is not 
closed downstream. At some distance from the breakdown point 
the flow again becomes more or less cylindrical. In other flows 
(on delta wings, for example) a stagnation point occasionally 
does not seem to exist; there is only a swelling of the core with 
an attendant decrease of velocities in the core. This situation 
is depicted in Fig. 1.4; the term 'vortex bursting* should 
probably be reserved for it. 

Finally, if the swirl is very high in a vortex, the vortex 
assumes a decidedly different columnar character where flow 
dependence on axial position has essentially disappeared, and 
the 'Taylor column' stretches far downstream (Fig. 1.5). This 
structure is found downstream of a vortex jump or immediately 
downstream of the flow entrance if the initial swirl is very 
high. In keeping with accepted notation (Benjamin 1962) this 
flow is termed ' subcritical ' while the flow upstream of a vortex 
jump would be 'supercritical'. (Stationary waves may occur in 
'subcritical' flow, while in 'supercritical' flow the axial flow 
is too fast, sweeping any waves downstream) . 

Summarizing the experimental observations, it is obvious 
that clear distinctions exist between 'vortex breakdown' and 
'vortex jump', 'supercritical' and 'subcritical ' flow, 'vortex 
breakdown* and 'vortex burst' , 'axisymmetric breakdown' and 
nonaxisymmetric instabilities. For convenience , we have also 
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Fief* 1*5 Subcritical (Taylor column) behavior. 
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distinguished between 'free' and 'forced' vortex bubbles, 
although they appear to be the result of the same basic, mecha- 
nism. In this report, the term 'vortex breakdown' will be used to 
designate conditions where the quasi- cylindrical solutions become 
singular and the corresponding computations fail. This includes rapid 
core expansion, as in the vicinity of the first stagnation point, and 
rapid core contraction. 

Witn accurate optical velocity measurements ahead of, 
around, behind, and in the vortex bubble now available (Orloff 
1971, Orloff and Bossel 1971) , it has become possible to test 
the accuracy and applicability of numerical approaches to the 
solution of the 'vortex breakdown' problem. It is found that 
the quasi-cylindrical approximation can be correctly applied 
upstream to a distance of the order of the bubble diameter ahead 
of the stagnation point (Figs. 1.1-1. 3). Flow in the neighbor- 
hood of the stagnation point and in and around the bubble must 
be treated by using either the full Navier-Stokes equations or 
a proper approximate subset (Figs. 1.1-1. 3). It was shown in 
Bossel (1967) that the proper subset is the equation for 
inviscid rotating flow. Numerical solutions show excellent 
quantitative agreement with the measured velocity profiles with 
an inaccurate description of the swirl velocity profile inside 
the bubble as the only exception (Orloff and Bossel 1971) . For 
accurate description of the very slow flow inside the bubble, 
viscous terms should be taken into account? however the exact 
bubble shape and the external flow are quite obviously deter- 
mined almost solely by the inviscid flowfield. The inviscid 
equation of rotating flow can be applied to compute 'free' and 
'forced' vortex bubbles equally well (Bossel 1967, 1969, Orloff 
1971, Orloff and Bossel 1971) . 

A continuous solution is not possible in a region where a 
vortex jump occurs; the 'flow force' concepts of Benjamin (1962) 
must then be applied to join two conjugate flows. In the 
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observations of Orloff and Bossel (1971) the vortex jump occurs 
over a distance of approximately one bubble diameter and is 
marked by very strong turbulence. 

1.2 Swirl Parameter 

The axisymmetric behavior of a vortex flow is very much a 
function of its swirl. A swirl parameter can be defined as 
the ratio of some representative (circumferential) swirl velocity 
to a representative axial velocity (see Sec. Ill) . For zero 
swirl, the flow reduces to a jet or wake, depending on whether 
a velocity excess or deficit exists on the axis. This jet or 
wake velocity profile will ultimately decay to the freestream 
velocity as the flow proceeds downstream. The same is true for 
small and medium amounts of swirl, until the swirl velocity and 
axial velocity reach comparable magnitudes . At higher swirl 
parameters, the asymptotic approach to freestream conditions 
is replaced by deceleration of the flow on the axis until a 
stagnation point and subsequent 'breakdown bubble 1 , or at least 
a swelling of the core appear. All other conditions being 
equal, there is obviously one particular swirl parameter value 
which divides a vortex flow which decays in a wake-like manner 
from another which 'breaks down'. This dividing swirl parameter 
is designated S 0 in the following; it is a function of the 
vortex velocity profiles. 

Axisymmetric vortex breakdown is a completely symmetric 
phenomenon with its origin at the very axis (Sarpkaya 1971a, b, 
Orloff 1971) . Thus conditions on the axis must be of vital 
importance in its inception. In particular, viscosity estab- 
lishes rigid rotation at and near the axis. For a rigidly 
rotating inviscid cylindrical slug of fluid of radius r , 
axial velocity u , and swirl velocity w at r , a swirl parame- 

dX C C 

ter can be defined by S = w /u , and two important theoretical 

c dx 

results follow: 
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(1) As the swirl is increased and crosses S Q = \/2 from 
below, an initial deceleration of the core element will 
amplify and lead to stagnation (Bossel 1968) . 

(2) As the swirl is increased further and crosses - hi'* 
(= 3.8317/2 = 1.9159)* from below, the character of 

the solution changes from 'supercritical' (which cannot 
support standing waves) to 1 Subcritical 1 (which can 
support standing waves) (Fraenkel 1956) . As S reaches 
the values j^ n /2 corresponding to higher zeros of the 
Bessel function J^, more waves appear. 

It can be expected that the critical swirl values S Q , S^, 

S 2 / . . . of a viscous vortex will be a function of the vortex 
velocity profiles and will be somewhat different from those for 
the rigidly rotating inviscid cylindrical slug. Of particular 
significance for vortex flows are S Q (dividing wake-type from 
breakdown behavior) and (beyond substantial upstream 
influence is possible) . 

1.3 Numerical Approach 

Laminar incompressible steady axisymmetrio vortex flows 
are described by the corresponding form of the Navier-Stokes 
equations. These equations can be approximated by a parabolic 
viscous set, analogous to the boundary layer equations, in 
regions where the stream surface angle remains small (quasi- 
cylindrical vortex flow) , and by the inviscid equations of 
rotating flow at and near the axis , and where stream surface 
angles become large (expansion or contraction of the core) 
(Bossel 1969) . 

Some solutions of the inviscid set pertaining to vortex 
flows at high swirl have been presented in Bossel (1967, 1969), 
Chow (1969) , Orloff and Bossel (1971) ; here we shall now discuss 

* 

is the first zero of the Bessel function J^. 
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a method of solving the parabolic viscous set and give corres- 
ponding results for different swirl parameters, initial velocity 
profiles, and external pressure and circulation gradients. 

Since the regions of validity of the two sets partially overlap 
on and near the axis , it should be possible to confirm some of 
the earlier results. The parabolic system has been used before 
in numerical computations by different methods (Gartshore 1963, 
Hall 1965, Bossel 1967, Mager 1971). Gartshore and Mager each 
used a momentum-integral approach and encountered a singularity 
which they linked to the vortex breakdown phenomenon and to the 
critical swirl parameter 3.8317/2 of flow in initially rigid 
rotation. The present results reinforce this view. Beyond 
this critical swirl ratio, Bossel (1967) and Mager (1971) also 
obtained flows which contrasted in behavior with those below 
thie critical swirl ratio. 

The numerical method (N-parameter 'exponential series 
integral method'*) to be used for the present study has 
previously been outlined in Bossel (1970a, 1971) . A related 
method (N-parameter 'power series integral method') has earlier 
been used successfully to calculate viscous vortex flows 
(Bossel 1967) . The methods were inspired by the success of 
the Dorodnitsyn method for the calculation of boundary layers 
(Dorodnitsyn 1962, Bethel 1968), especially by its aspects of 
speed and accuracy. However, the Dorodnitsyn approach with its 
use of Y (U) in place of U(Y) cannot be directly applied to 
vortex flows where axial velocity profiles often exhibit over- 
shoot and flow reversal. The present approach circumvents this 
problem. The soundness of the exponential series integral 
method has been demonstrated in its application to incompressi- 
ble and compressible boundary layer flows (Bossel 1970a, b, 

Mitra and Bossel 1971) . 

* 

The term 'integral method' is used here to avoid the cumber- 
some (but more accurate) terms 'method of weighted residuals' 
or 'method of integral relations'. 
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There are good reasons for choosing an N-parameter integral 
method over a finite difference method for a study such as the 
present one. Many of the open questions in this case only 
require a qualitative answer, which can be obtained at minimal 
computing cost by a one or two parameter solution. If accuracy 
is desired, the number of parameters is increased without any 
change in the program. Results of the same detail and accuracy 
as for finite difference methods can then be obtained. The 
integral approach reduces the solution of partial differential 
equations to solution of a set of ordinary differential equa- 
tions. Well-known stable integration schemes, such as the 
Runge-Kutta method, can be applied. There is little chance of 
confusing a true singularity of the set of equations to be 
solved with a singularity of a finite difference scheme replacing 
it. No iterations are necessary. Finally, velocities, stream 
function, pressure, boundary layer thickness, shear, etc., can 
all be computed simply and directly from analytical expressions 
involving the parameters in the velocity approximations. 

1 . 4 Overview 

The following Section II will first describe the computa- 
tional method and give an example for its accuracy and 
convergence. Singularities of the equations are investigated in 
Section III, and a swirl parameter is introduced. Section IV 
investigates , for constant external axial velocity and circula- 
tion, three vortex flows which are of particular interest to 
the fluid dynamicist: (1) vortex flow with initial uniform axial 

velocity, (2) vortex flow where the velocity on the axis is 
initially higher than the frees tream velocity (leading edge 
vortex) , and (3) vortex flow where the velocity on the axis is 
initially lower than the freestream velocity (trailing vortex) . 
Representative velocity profiles and velocity distributions on 
the axis as a function of distance are given. The effects of 
positive or negative gradients in the external velocity and 
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circulation distribution are studied in Section V. Some 
conclusions are drawn, and results are summarized in Section VI. 
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II. COMPUTATIONAL METHOD 


In developing a method of weighted residuals for a problem 
in two coordinates, the dependent variables in the partial 
differential equations are replaced by approximating functions, 
multiplied by a set of weighting functions, and then formally 
integrated with respect to one coordinate. There remains a set 
of ordinary differential equations in the other coordinate for 
the parameters in the approximating expressions for the dependent 
variables. The major task in developing the integral method is 
the correct development of the coefficients in the set of ordinary 
differential equations - not a difficult task, but one requiring 
careful bookkeeping. The FORMAC* method has been used in a case 
of comparable complexity to relegate this task to the computer 
(Mitra 1970) . The ordinary differential equations are solved for 
the parameters by any of a number of standard methods. Knowledge 
of the parameters at each step permits calculation of the depen- 
dent variables and of a number of related variables of interest. 


2 . 1 Governing Equations 

The independent and dependent variables are introduced in 
Fig. 2.1. The dimensional equations for incompressible viscous 
quasi-cylindrical vortex flow are (Bossel 1969) : 


3u ■ 3v , v 
3x 3r r 


0 


W 2 _ 1_ 3£ 
r p 3r 


u 


u 


3w 

3x 



3u 

3r 

i |E + V 
p 3x 

/3 2 u 

Ur 2 

+ i 

r 

3u^ 

3r/ 


vw 

— " v 

1 + I 

3w 



L 2 r 

3r 

r 2 / 



( 2 . 1 ) 


* 

FORmula MAnipulation on Computer , a method developed by IBM 
for non-numerical manipulations. It uses PL-1 or FORTRAN. 
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Fig. 2.1 Coordinate system and velocities. 


with the boundary conditions: 

at r = 0 : v = 0, w = 0 , u = u (free) 

at r + »: u -*■ u (x) , wr = k k (x) 

e e 

Nondimensional variables are defined as follows : 

X = x/r c Y = R 2 /2 = Re(r 2 /2r c 2 ) P = p/(pu oo 2 /2) 

U = u/Uoq V = /Re (v/u OT ) W = w/u^ 

where the core Reynolds number Re = u^r /v . For convenience we 
define H = VR and K = WR (circulation) . These quantities are 
introduced into the conservation equations. The pressure P is 
eliminated by cross-differentiation, and two equations remain 
for U and K : 
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w <UK> + A ( hk - 2y # + 2k ) - ° 

W ( 2y2u w + t) + y2 ^T ( hu - 2Y W ) - 0 

where Y 

0 

The boundary conditions now become: 

at y = 0: H = 0 , K = 0 , U=U (free) 

at Y -*■ * : U -> U g (X) , ?7R = K ■> K & (X) 


The pressure at any point in the flow is given by 


P = P 

e 



dY 


( 2 . 2 ) 


where 



U 

e 


2 


2 . 2 Integral Relations 

As a first step in developing the computational method, 
integral relations are derived from the two momentum equations 
(2.2). In order to obtain a sufficient number of equations for 
determination of the unknown parameters in the axial velocity 
and circulation approximations to be introduced below, the 
equations are multiplied by members of two sets of weighting 
functions g^(Y) and f^ (Y) , respectively, and then integrated 
formally in the Y-direction from zero to infinity. To facilitate 
the integration, we require g fc (0), ( 0 ) = finite, and g k (”) , 

f k (°°) = 0. After integration by parts, rearrangement, and 
simplification : 
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dX 


00 oo oo 

J g k UK dY - f g k ' HK dY - f (2g k "Y + 4g k ')K dY = 0 
0 0 0 

00 °° 2 

4- J (f k 'Y 2 + 2f k Y)U 2 dY - / f k dY 

0 0 

00 

- j (f K "Y 2 + 4Yf k + 2 f k )HU dY 
0 

00 

- f ( 2 f k Y 3 + 14f k "Y 2 + 2 0 f k ' Y + 4f k )U dY = 0 

0 


(2.3) 


2 . 3 Approximating and Weighting Functions 

The integration is completed with the following choices for 
weighting functions and axial velocity and circulation approxi- 
mations : 
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• 



The external velocity and circulation U (X) and K (X) are 

w G 

prescribed; the a (X), b (X) and the velocity on the axis U (X) 

n n ax 

are free parameters. Major reasons for these choices (Bossel 
1970a) are the exponential character of the flows, ease of 
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analytical integration , the fact that the approximations satisfy 
the boundary conditions , and the observation that the approxi- 
mations satisfy Weierstrass' approximation theorem after a 

coordinate transformation of the semi-infinite region 0 Y < 00 

— rv Y 

into 1 > n > 0, with n - e 


2.4 System of Ordinary Differential Equations 

Introduction of the approximating expressions and the 
weighting functions (2.4) into equation (2.3) results in the 
following set of ordinary differential equations for U „ and the 

3.X 

parameters a (X) and b (X): 
n n 


N N 

2 a A . + > b B . + U C. = -U 

n n,k £_> n n,k ax k 

n=l 


D, - K E, - F, 
e k e k k 


n=l 


k = 1, 2, . . . , N + 1 

2 A nVk + 2 Vn,k + 6 aA “ -°e 5 k ' KK * f k 

k = 1 , 2, . . . , N 


N 


n=l 


N 


n=l 


(2.5) 


The coefficients are given in the Appendix. Only a small part 

of the coefficient calculation has to be done at each step; the 

major part is done only once at the beginning of the computation 

The system (2.5) of (2N + 1) first order ordinary differential 

equations is first solved (by Gaussian elimination) for the 

derivative vector (da /d X; db /dX; du /dx) . A standard method 

n n dx 

of integration then produces the a n , b R , and U ax (the Runge- 
Kutta method has been used in the present work) . Axial velocity 
U(X f Y) , circulation K(X,Y) , and swirl velocity W(X,Y) follow 
from relations (2.4). 
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2.5 Initial Profiles 


Initial profile parameters a , b , and U , for the initial 

n n ax 

axial velocity and circulation profiles (2.4) are required to 
start the calculation. In the present study, simple exponential 
profiles have been used throughout, i.e.. 


“initial' tY > ~ V 1 - e '“ Y > + a ax e '“ Y 


and 


“initial (Y > - 


K (1 - e -otY ) or K (1— e" 2olY ) 
e e 


Exponents a = 1, and = k, k = 1, 2, . . ,,N+1 were used in 

all of the calculations. 


In the case of arbitrary profiles, the parameters a n , b n , 
and U are obtained as outlined in Bossel (1970a or b) . The 

aX 

procedures for the axial velocity profile U(Y) and the circula- 
tion profile K(Y) are identical. Thus, if U(Y) is the 
(analytical or tabulated) profile to be approximated, U N (Y;a ) 
the profile approximation (with unknown parameters a n ) , and 
f^(Y) a convenient weighting function, then the N unknown 
parameters a^ are obtained from the N equations 


/ f k (Y)[U(Y) - U N (Y;a n )] dY = 0 


N 


2 . 6 Accuracy and Convergence 

Results of the present method were compared with previous 
computations (Bossel 1967) using polynomials in Y in the approxi 
mating functions. For N = 2 the agreement was within 3% in the 
velocity profiles. In cases where they can be compared, the 
computations appear to agree well with those of Hall (1966b) 
(only a qualitative comparison can be made, however, since Hall 
prescribes the outer stream surface shape) . Critical runs have 
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been repeated with different orders N of approximation. Figures 
2.2 and 2.3 show the velocity on the axis, swirl parameter, and 
profile development for a breakdown case (type lb) and N = 1, 

2 , 3 , 4 , and 5 . There is obvious convergence as N increases . 
Cases like these have typical running times on the IBM 360/75 of 
5 to 15 seconds for N = 1 and 2, and 10 to 30 seconds for N = 3 
to 5. All cases discussed later were run with N = 3, which 
appeared to be an efficient compromise between the conflicting 
demands of high accuracy and low computing time. 



Fig. 2.2 Development of velocity on the axis and local swirl 
parameter S (X) as a function of number of parameters 
used, for initial swirl parameter S. = 1.425. 
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III. SINGULARITIES AND CRITICAL SWIRL PARAMETERS 

3 . 1 Singularities of the Governing System 

The system of ordinary differential equations (2.5) becomes 

singular when the determinant of coefficients for the unknowns 

a , b , and U vanishes , and the right-hand side remains 
n n cix 

nonzero. The special case when the right-hand side vanishes 
also has some significance and is considered later. 

At a particular axial position X Q the profile shapes 
(relations 2.4) are determined by the local a n (X Q ), b n (X Q ) , and 
U (X ) . Without loss in generality we shall assume U = 1 and 
keep only K as a free parameter. This amounts to investigating 
a given axial velocity profile and a given circulation profile 
shape at different values of external circulation K e> The 
condition for vanishing of the determinant of coefficients on 
the left-hand side of system (2.5) is represented by a polynomial 
expression of the form 

N n 

X G n (K e 2) = 0 (3.1) 

n=0 

with N roots for characteristic external circulation values 
2 

(K ). The G„ are numerical coefficients. Thus, there are N 
e n 2 

distinct values of K which lead to singularities, and for 
N ■+■ «> the following results obtains 

For a given axial velocity profile and given shape of the 
circulation profile a the equations of quasi- cylindrical viscous 
incompressible vortex flow (2.1) have a countably infinite set 
of singularities for discrete values of external circulation , 



It should be stressed that the singularities are a function 
of the local profiles only and independent of local axial gra- 
dients of external axial velocity and circulation. However, the 
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profile development itself is, of course, dependent on external 
gradients . 

The effect of gradients of external axial velocity U e and 
circulation K g (i.e. dU e /dX and dK e /dX) becomes obvious if the 
special case is considered where the right-hand side of equations 
(2.5) vanishes together with the determinant of coefficients . 
System (2.5) is then nonsingular. The right-hand side is a 
linear function of the gradients dU e /dX and dK e /dX and the 
parameters a , b , and U in the axial velocity and circulation 
profiles. For each profile combination (i.e. U g , K , a n , b , 

U ) a particular combination of external axial and circulation 
gradients dU /dX and dK /dX exists for which the right-hand 
side vanishes and the system (2.5) is nonsingular, i.e.: 

A singularity can be avoided by application of appropriate 
external axial and circulation gradients . 

This result is supported by the results of computations 
for different external axial and circulation gradients described 
in this report. 

The two results derived above are independent of the 
particular approximation and weighting functions chosen (i.e. 
expressions (2.4)) since any continuous function U (Y) or K(Y) 
can be represented by an infinite series of linearly independent 
functions {4> n (Y)} forming a complete set. 

It becomes evident from the analysis that the (swirl- 
dependent) singularities of the quasi-cylindrical vortex equations 
are different in character from the separation singularity of the 
boundary layer equations. In particular, in contrast to boundary 
layer solutions , there exist meaningful vortex solutions for 
swirl values between the singular values, and especially for a 
swirl value greater than that of the first singularity. Full 
vortex computations (Secs. IV and V) have shown, however, that 
these solutions will normally only persist for a finite distance, 
before they are driven to the next singularity. 
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3.2 Swirl Parameters 


In comparing different vortex flows, K g is not a very 
meaningful parameter. Since the axisymmetric breakdown behavior 
of vortex flows is evidently controlled by conditions at the 
very core (Sarpkaya 1971a, b, Orloff and Bossel 1971) , a swirl 
parameter is introduced which is based on core conditions, i.e. 
the swirl velocity profile gradient on the axis, the axial 
velocity on the axis , and the core radius R where the 

swirl velocity W reaches its maximum value: 


S = 


(dW/dR) 

U 


ax w max 


ax 


This definition of the swirl parameter reduces to the familiar 

one S = ftR /U = W /U for a cylinder of fluid of radius R in 
c c c 

rigid rotation ft , with uniform axial velocity U. Inviscid 
solutions for such flows (Fraenkel 1956) have critical values 
of the swirl parameter at j^ n /2 = 1.9159, 3.5078, 5.0867, 
6.6618, 8.2353, ..., where j^ n are zeros of the Bessel function 


3 . 3 Singularities for Initially Uniform Axial Flow 

Figure 3.1 presents the first five singularities for the 
vortex profiles 

U (Y) = 1 = constant 

K ( Y) = WR = K (1 - e“ aY ) = S (1 - e -aY )/0.792 

as a function of swirl parameter and order of approximation used. 
This flow approximates rigid rotation in the core. The singular- 
ities were obtained by computing the derivatives a n , £> n , U ax 
over a wide range of K e> Each additional approximating term 
accounts for a new singularity, as expected, while the location 
of previously computed singularities is only slightly affected. 
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0 1 2 3 4 5 inviscid 

N rigid rotation 

Pig. 3.1 Singular swirl parameter values for U (Y) = 1, and 
K(Y) = WR = S (l-e” Y ) /0 . 792 as a function of order 
of approximation N. (Inviscid values are for rigid 
body rotation.) 



Fig. 3.2 Magnitude of gradient of velocity on the axis as 
function of swirl parameter. Uniform axial flow 
U (Y) = 1, and K (Y) = WR = S (1-e-Y) /0 . 792 . 




Likewise, the use of different a or weighting function exponents 
k has only minimal effect on the location of the lower order 
singularities. In particular, the singularity of greatest 
interest, S^, remains practically constant. Since there is 
ample evidence that this first singularity is associated with 
axisyipmetric vortex breakdown, computations failing at singular- 
ities, and in particular at S^, will also be said to 'breakdown' 
in the following work. The critical swirl values j^ n /2 as 
found from the theory of inviscid rigid rotation (Fraenkel 1956) 

' are also shown in Fig. 3.1. This figure also gives the swirl 
parameter S Q dividing flow which decays smoothly (type la) 
from flow which is abruptly decelerated on the axis and 'breaks 
down' because is reached (type lb) . The theoretical 
(inviscid) value (Bossel 1968) for uniform flow in initially 

rigid rotation is S = /2. The viscous value for is found 

o o 

from a full vortex computation (see Sec. IV) . There is striking 
agreement between the viscous and inviscid values for S Q , S^, 

S 2 , ... . 

\ 

3 . 4 Behavior of Axial Derivatives 

At each singularity, all axial derivatives jump to infinity 
in magnitude and reverse signs. This results in contrasting 
behavior for flows separated by a singularity. The axial 
derivative of the velocity on the axis for the vortex with uni- 
form axial flow is plotted in Fig. 3.2 as a function of swirl 
parameter S. The singularities are very evident. Note especial- 
ly that the sign of the axial gradient reverses at each 
singularity. In particular, for a swirl value below S^, the 
axial flow is decelerated in the core, while for a swirl value 
S-^ < S < S 2 the axial flow is accelerated in the inner core. 

3 . 5 E ffect of Velocity Profile 

The critical swirl parameters S Q , S^, S 2 , . are functions 
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of the a , b , U , and U and are therefore profile-dependent. 
The example given, which approximates rigid rotation in its 
inner core, appears to confirm the earlier contention (Bossel 
1967, 1968) that the behavior of viscous vortex flows is governed 
mainly by inviscid core properties. The effect of nonuniform 
initial axial flow on the critical swirl parameters is shown 
in Fig. 3.3 for the family of profiles 

U (Y) = (1 - e" aY ) + U e" Y = 1 + e” aY (U - 1) 

dA aA 

K (Y) = WR = K (1 - e“ aY ) 

Flows with U <1 are of the trailing vortex type, with a 

clX 

velocity deficit on the axis, while flows with U > 1 approxi- 
mate the leading edge vortex with a velocity excess on the axis. 

Figure 3.3 shows the same qualitative behavior for the 
various flows. The magnitudes of the critical swirl values for 



Fig. 3.3 Effect of velocity profile on the critical swirl 

values. U(Y) = 1 + e“ Y (U ax -1) , K (Y) = S (l-e" Y ) /0 . 792 . 
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the nonuniform profiles show a profile effect: they are 

consistently higher for the profile with initial axial velocity 
deficit, and consistently lower for the profile with axial excess 
near the axis. This result is to he expected. The average velocity near 
the axis is most certainly higher than U in the trailing vortex case 

£LX 

due to the effect of faster flow in a neighborhood of the axis. In the case 
of the leading edge vortex, fluid surrounding the axis is slower than 
and the average velocity near the axis would be less. If proper 

average velocities were used in the calculation of S , the lines of 

critical swirl parameters would be more nearly horizontal. 

Singularity analyses can be performed for any given arbi- 
trary (experimental) axial and swirl velocity profile combination 
in order to obtain information on likely flow behavior, in 
particular the proximity of S^. The coefficients a n , b n , and 
U in an N-th order profile approximation (2.4) are first 
determined by the method of Bossel (1970a) . The critical swirl 
values are then either obtained by determining the coefficients 
G n in the characteristic equation (3.1) and solving for the 
first N critical (or S cr ) , or by determining the gradients 
d/dx of the a n and b n from system (2.5) as a function of K g (or 
S ) . 


3 . 6 Assessment of Breakdown Behavior 

The results of this section offer several possibilities 
of assessing the probable axisymmetric breakdown behavior of 
vortex velocity profiles. In the order of increasing accuracy: 

(1) Comparison of the swirl value 

„ . <dW/dE) ax R c _ <dw/dr) ax r c 
D ax 

with the inviscid critical values for rigid rotation, in 
particular S Q = fl (dividing stagnating (S > S Q ) from 
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wake-type vortex flow (S < S Q ) ) , and = 1.9159 
(dividing supercritical (S < S^) from subcritical 
vortex flow (S > S-^) , with a reversal of behavior at 

. 

(2) Comparison of the computed swirl value with the results 
of Fig. 3.3. 

(3) Approximation of the velocity profiles by parameters 

a , b , U , U , K using the method in Bossel (1970a) 
n n e ax e 

for initial profile approximation. Using the approach 
of the present section, the axial gradients are then 
computed and the actual singularities determined by 
varying K g or the swirl S . 
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IV. VORTEX COMPUTATIONS FOR CONSTANT EXTERNAL AXIAL VELOCITY 
AND CIRCULATION 


This section presents results of vortex computations for 
constant external axial velocity and circulation. Three 
different initial velocity profiles were used in this and in the 
following section, corresponding to initially uniform axial flow 
and to flows of leading edge and trailing vortex type, respec- 
tively. Diagrams of distance X^ to failure of the computation 
as function of initial swirl parameter will first be presented 
for these three flow types. The development of the velocity on 
the axis and of the velocity profiles will be discussed in some 
detail. The results confirm the singular behavior predicted in 
Sec. III. 

4 . 1 Initial Profiles and Distance to Failure 
Uni_foj^rn ^initial axial 

The case where the initial axial velocity profile is uniform 
while the swirl velocity profile is linear with R in the vicinity 
of the axis provides a good test of the assertion that vortex 
behavior is governed by the inviscid properties of the rigidly 
rotating core as conditions on the axis are then identical in 
both cases. Initial profiles (Burgers 1940) for this case were 

U(Y) = 1 = constant 

K (Y) = W (Y) R = K (1 - e -yY ) 

Most computations presented here were obtained for y = 1. 
y = 2 was used to test scaling effects on the results of the 
computation. The differences were found to be minimal. The rela 1 
tionship between external circulation and initial swirl can be 
computed from the initial profiles. For y = 1 
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I I 


s = 


<aw/dR) ax R c 


u 


ax 


- - 792 V° ax 


.792 K 

e 


and for y = 2 

S = 1.121 K /U = 1.121 K 
e ax e 

Figure 4.1 presents, as a function of initial swirl para- 
meter and for N = 3, the distance of successful computation 
before failure occurred at X^. This plot shows several distinct 
regions separated by singular points in accordance with the 
results of Sec. III. Continuous solutions could only be found 
in region la, corresponding to an initial swirl parameter 
S 5 /2. The characters of the solutions in the different regions 
will be discussed more fully in the next section. As pointed 
out in Sec. Ill, the computation with N = 3 only picks up the 
first three singularities. The results are therefore only 
accurate for S $ S 2 . 



Fig. 4.1 Distance to failure (at Xf) as function of initial 
swirl parameter. Initially uniform axial flow. 
Open circles denote failure of computation. 
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Leading edge vortecc 


Experimental (Hummel 1965) and theoretical (Hall 1961) 
results for leading edge vortex flows (as generated at the 
leading edges of delta wings, for example) all show a large 
excess of axial velocity near the axis over the freestream axial 
velocity. It is of some interest to see what effect, if any, 
this velocity excess has on the vortex behavior. The effect of 
the nonuniform shape of the axial velocity profile on the 
critical swirl values has been discussed in Sec. III. 

The initial velocity profiles investigated were 

U (Y) = 1 + e“ Y (i.e. U = 2.0, U = 1.0) 

aX 6 

K (Y) = W(Y)R = K (1 - e“ Y ) 

The relationship between K g and S for these profiles is 

S = .792 K /U = .792 K /2 = .396 K 
e ax e e 

This choice is not too good an approximation to the experimental 
and theoretical profiles of Hummel and Hall, but it should serve 
to demonstrate the effect of higher axial core velocity. 

Distances to failure of the computation at for N = 3 
are shown in Fig. 4.2 as a function of initial swirl parameter. 
The behavior is qualitatively the same as for the case of initial 
uniform axial flow (Figure 4.1). The critical swirl values are 
now somewhat different from the previous case. They have been 
more accurately determined, with N = 5, in Sec. III. 

Trailing vortex 

The trailing vortex, as found downstream of a straight 
wing, has a wake-like axial velocity distribution with a core 
velocity smaller than the freestream velocity (Batchelor 1964, 
McCormick et al . 1968). This velocity retardation near the 
axis can be expected to lead to earlier vortex stagnation and 
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, „ breakdown- stable 

0 I I I I l 1 — l 1 — „.i . I , J 
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X f 

Fig. 4.2 Distance to failure (at Xf) as function of initial 
swirl parameter. Initial profiles of leading edge 
vortex type. Open circles denote failure of 
computation . 


breakdown (see Sec. Ill) . 

The initial profile chosen was 


U (Y) = 1 - e /2 


(i.e. U. 


0.5, U e = 1.0) 


K (Y) = K e (l - e x ) 


For these profiles, external circulation and swirl values are 
related by 

S = .792 K /U = .792 K /0.5 = 1.584 K 
e ax e e 


This choice appears to come fairly close to actual trailing 
vortex flows (McCormick et al. 1968) . 

Distances to failure of the computation at X^ for N = 3 are 
shown in Fig. 4.3 as a function of the initial swirl parameter. 
Again the behavior of the solutions is qualitatively the same as 
for the two previous cases. The actual critical swirl parameter 
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Distance to failure (at Xf) as function of initial 
swirl parameter. Initial profiles of trailing vortex 
type. Open circles denote failure of computation. 


values differ somewhat from the other two cases and again are 
not as accurate as those found in Sec. Ill for N = 5, especially 
for S > S£. 

o£ breakdowri^ solutions of t}j] 2 _e 2 3 g 3 4^ ... 

No breakdown-stable solutions of type 2, 3, or 4 were 
found in the present investigation for constant external axial 
velocity and circulation. This agrees with the results of the 
singularity analysis in Sec. Ill where it became evident that 
only specific external axial velocity and circulation gradients 
can prevent flows of these types from approaching a singularity. 
It will again be shown in Sec. V in several examples that 
application of external axial velocity or circulation gradients 
can have beneficial effects on flow development for these types. 
There is some other experimental and analytical evidence 
(Donaldson and Sullivan 1960) suggesting breakdown-stable solu- 
tions of the multi-layered type. 
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4.2 Velocity on the Axis 

The behavior of vortex flow is different for each of the 
regions la, lb, 2, 3, ... . The most important indicator of 
vortex behavior is a plot of velocity on the axis U vs. 
distance X. The development of the velocity on the axis will 
now be discussed for flows with constant external axial velocity 
and circulation. Before describing in greater detail flows of 
type 1 having different initial profiles, behavior of initially 
uniform axial flows of the first five types will be discussed. 
All results were computed with N = 3. 

B_ehavioZ. °f axial flows of t.ypes 1-4 

Figure 4.4 presents the development of velocity on the axis 
in the five regions la, lb, 2, 3, and 4. Each of the curves is 
representative of other flows in the same region. For flows of 
a given type, the character of these curves, and of the velocity 
profiles, remains the same. Type la, the only breakdown-stable 



Fig. 4.4 Development of velocity on the axis for different 
flow types. Initially uniform axial flow. Open 
circles denote failure points. 
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flow, develops a reduce velocity on the axis which remains 
approximately constant for a considerable distance before it 
again approaches the freestream value. Type lb flow has 
gradients (dU /dx) which are always steeper than any type la 
gradients. These gradients increase, and the computation 
terminates due to excessive gradients, with U in the neighbor- 

3.X 

hood of 0.5. At the same time, the stream surface angles 
suddenly increase to large positive values, indicating an 
explosive behavior of the core. Gradients become larger, and 
the calculation breaks down increasingly sooner, as the critical 
swirl value S^ is approached. 

For swirl values greater than S-^, solutions of type 2 are 
obtained, which now show an acceleration of the velocity on the 
axis. Gradients (dU /dX) are very large for S near S. , but 
decrease as S is increased, leading to longer calculations before 
failure occurs. Failure of type 2 is an implosive one: stream 

surface angles assume large negative values very quickly in the 
immediate vicinity of the failure point. Gradients decrease 
further with correspondingly longer resulting computations, 
until the vicinity of S 2 is reached, where again no solutions 
can be obtained. 

For swirl S > S 2 / the flow is again of different type 3. 

The velocity on the axis now again decelerates , eventually 
reaching a stagnation point on the axis. Downstream of this 
point reversed flow develops on and near the axis, gradients 
steepen and computation fails with explosive behavior of the 
outer core. As the swirl is further increased, gradients are 
again reduced, and the computations fail farther downstream 
until the vicinity of S 3 is reached. No solutions can be ob- 
tained in the immediate neighborhood of this point. 

For swirl S > S 3 , the flow has once again changed its 
behavior (to type 4 flow) . The velocity on the axis again tends 
to increase, until gradients again steepen rapidly, leading to 
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large negative stream surface angles in the outer core regions 
(implosive behavior) , and the computation fails. 

As S is increased, more and more singularities are encoun- 
tered. The velocity on the axis alternates between a decelerating 
tendency (for odd type flow) and an accelerating tendency (for 
even type flow) . 

Initial^ ani/orm axial /Zoft? of tijge 1 

Probably the only vortex flow of major physical signifi- 
cance is type 1 flow. The development of its velocity on the 
axis and its swirl parameter S(X) will now be given in more 
detail for the vortex with initially uniform axial flow. 


Figure 4.5 presents results for different choices of initial 
swirl parameter S^. It also gives the swirl parameter S (X) for 
two neighboring flows of type la and lb. A swirl parameter 


value of approximately = 1.4 separates the smoothly decaying 
vortex flows (type la: £ 1.4) from vortex flows which 


Fig. 



4.5 Development of velocity on the axis as a function of 
initial swirl. Initially uniform axial flow. Open 
circles denote failure of the computation. 


36 



eventually break down (type lb: £ 1.4). Note that the 

theoretical inviscid value for the initial swirl value leading 
to stagnation is > /2 = 1.41 for rigid rotation (Bossel 1968). 

Type la flow adjusts to an almost constant value of axial 
velocity and thereafter behaves in a wake-like manner. The 
swirl parameter S (X) generally decreases. The swirl still 
counteracts the normal spreading of the wake and the attendant 
diminishing of the axial velocity deficit? but eventually the 
axial velocity deficit again decreases , and the axial and swirl 
velocity profiles decay together. 

Flows of type lb initially behave much like those of type 
la, except that their rate of decrease of U is steeper, and 
the swirl parameter S (X) increases. Eventually the swirl effects 
overwhelm the restoring tendencies of the wake, the drop of 
velocity on the axis steepens rapidly, and the computation 
fails. As an initial swirl parameter value of a 1.8 is 
approached, the computation fails increasingly sooner. No 
solutions can be obtained near this point, whose theoretical 
inviscid value is = 3.8317/2 = 1.9159 for rigid rotation 
(Fraenkel 1956) . 

Leadiri£_ edge cortex of_ tj£2S 1 

The development of the velocity on the axis for different 

initial swirl parameters for a vortex of leading edge type 

(initially U = 2, U = 1) is given in Fig. 4.6. The swirl 

parameter S(X) is also shown for two neighboring flows of type 

la and lb. The behavior is qualitatively again the same as for 

the vortex with initially uniform axial velocity. A swirl 

parameter value of = 1.1 divides the smoothly decaying type 

la from the breakdown type lb. This value of S. is somewhat 

lower than the theoretical value S = for inviscid flow in 

o 

rigid rotation and appears to reflect the influence of the 
profile shape. One would expect that the regions adjacent to 
the axis, where the axial velocity is lower, would have some 
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Fig. 4.6 Development of velocity on the axis as a function of 
initial swirl. Initial flow of leading edge vortex 
type. Open circles denote failure of the computations. 


effect on the flow behavior and would therefore tend to bring 
down the critical value of S^. 

TraiZ-inff o_f_ tyj^e 1 

The development of velocity on the axis for different 
initial swirl parameters Si for a trailing vortex (initially 
U = 0.5, U = 1) is plotted in Fig. 4.7. The swirl parameter 
S(X) is also given for two neighboring flows of type la and lb. 
The dividing value S q of the initial swirl parameter is now 
greater than /2 , reflecting again the influence of the nonuni- 
form initial axial velocity profile. 

4 . 3 Velocity Profiles 

Figure 4.8 presents typical velocity profiles for initially 
uniform axial flow and N = 3 as they develop in the different 
swirl regions. The swirl velocity profile is hardly affected, 
but the axial velocity profile shows characteristic features 
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X 


Fig. 4.7 Development of velocity on the axis as a function 
of initial swirl. Initial flow of trailing vortex 
type. Open circle denotes failure point. 


Type 1 Type 2 Type 3 Type 4 

Sj =1.39 Sj =3.36 Si =7.85 Sj =16.81 

X=0.04 X=0.06 X =0.074 X=0.08 



Idealized 

stream surface pattern 
near failure 




Fig. 4.8 Velocity profile development and flow character for 
different flow types. Initially uniform axial flow 




which are different in each region. Axial velocity profiles for 
these types, and for the initial profiles with uniform axial 
velocity, and leading edge vortex and trailing vortex character 
are given in more detail below. 

For initially uniform axial flow a type 1 vortex always 
shows a decrease of velocity on the axis while in type 2 flow, 
the increase in velocity on the axis appears, as noted earlier. 
In addition, flow in a layer adjacent to the core may now be 
retarded, and the vortex may assume a two-layer structure. 
Similarly, in type 3 flow, three layers appear: a retarded 

core flow, an adjacent accelerated layer, and again, an outer 
retarded layer. In type 4 flow, four similar layers appear, 
with core flow being accelerated. This core structure has been 
qualitatively confirmed by calculations for different numbers 
of parameters N in the integral method of solution. 

The velocity profile development for a smoothly decaying 
type la vortex with initially uniform axial flow is shown in 
Fig. 4.9. Behavior of type lb flow is qualitatively similar up 
to the failure point. The swirl velocity profile shows a 
gradual decay, while the effect on the axial velocity profile 
(axial velocity retardation on and near the axis) is more pro- 
nounced and increases with swirl. Figures 4.10, 4.11, and 4.12 
show the development of velocity profiles (y = 2) when the 
initial swirl parameter exceeds S^. These three cases of 
types 2, 3, and 4 are representative of families of solutions 
between singularities and S^. These solutions develop a 
multi-layer structure, and the core flow reacts differently in 
each family (i.e. acceleration or deceleration on the axis). 

All of these cases eventually fail. 

Figure 4.13 shows the development of the axial and swirl 
velocity profiles for a (breakdown) type lb flow of leading edge 
vortex type. Again there is a strong retarding effect on the 
axial velocity on and near the axis, while the swirl velocity 
profile is not greatly affected. 
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Fig. 4.12 


Velocity profile dev 
axial flow. Type 4. 


for initially uniform 







Fig. 4.13 Velocity profile development for initial flow of 
leading edge vortex type . Type lb . 


Figure 4.14 presents the development of axial and swirl 
velocity profiles for smoothly decaying type la flow of trailing 


vortex type . 


S| = 1.864 


Velocity Scale 



Fig. 4.14 Velocity profile development for initial flow of 
trailing vortex type. Type la. 




V. EFFECTS OF EXTERNAL VELOCITY AND CIRCULATION DISTRIBUTION 

The strong effects of external axial pressure gradients 
on vortex flows are well established experimentally (Lambourne 
and Bryer 1961, Harvey 1962, Hall 1966a) and in numerical 
computations (Hall 1966b, Mager 1971) . A positive external 
velocity gradient has an accelerating effect on the axial velo- 
city on and near the axis in type 1 vortex flows , thus decreasing 
the local swirl parameter and delaying or preventing breakdown. 
The effects of circulation gradient are less obvious. Combined 
external axial velocity and circulation effects may reinforce or 
cancel each other, depending on sign and magnitude of each and 
on the type of vortex flow. It will be shown here that external 
axial velocity and circulation gradients have contrasting effects 
on flows of different type. Such gradients can be used to 
prevent vortex breakdown, at least for some distance. 

The behavior of vortex flows under external velocity and 
circulation gradients is of considerable interest in cases 
where breakdown of vortex flows is to be prevented or induced. 
Computation by methods such as the present one can indicate 
what velocity and circulation gradients should be applied and 
where . 

5 . 1 Type 1 Vortex Flow 

Type 1 flow is of considerable interest because it is the 
most common type in applications. Leading edge and trailing 
vortices on wings are usually of the breakdown-stable type la, 
but may change to stagnating type lb as a result of either 
higher swirl (increased wing angle of attack) or adverse pressure 
gradient (negative external axial velocity and/or external 
circulation gradient ) , 

External axial velocity and circulation gradients are felt 
through their effects on both the external pressure distribution 
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and the* swirl parameter . The Bernoulli equation for the free- 
stream 

P Q - P e = (P/2) (u e 2 + v e 2 + (k e / r ) 2 ) 

« (p/2)Cu e 2 + (k e /r) 2 ) 

shows the same qualitative effect of external axial velocity and 
circulation gradients on the pressure gradient. There is an 
amplified response of the flow on and near the axis to external 
pressure gradient (see Hall 1966a, p. 69) which is confirmed by. 
all the present computations. Thus positive external axial 
velocity and circulation gradients will partly tend to cause 
increasing velocity on the axis for type 1 flow. 

Flow behavior is governed by the swirl parameter S ~ K /U 
While the effect of positive external axial velocity and 
circulation gradients increases U ax and decreases S, a circula- 
tion gradient also increases K and may raise S. A positive 
external velocity gradient is therefore always stabilizing in 
type 1 flow, while a positive external circulation gradient 
may be either stabilizing (in type la flow) , or destabilizing 
(in type lb flow), depending on which effect dominates.* The 
results presented here are for type 1 vortex flows with initially 
uniform axial velocity. Other initial axial velocity distribu- 
tions give qualitatively the same result. 

Figure 5.1 presents the results of the application of 
various positive and negative external axial velocity gradients 
to a type lb vortex flow of fixed initial swirl and with 
initially uniform axial velocity. This particular flow breaks 
down under zero external velocity gradient. The same qualitative 
behavior is found for a type la vortex flow. The experimental 

* 

'Stabilization' is here taken to mean the prevention of axi- 
symmetric vortex breakdown by avoiding any of the 
singularities S^, S 2 r ••• 
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Fig. 5.1 



of external axial velocity gradient. Initially 
uniform axial flow of type lb with = 1.425 and 
K e = 1.8 = const. 


observations (breakdown delay or avoidance through a positive 
external axial velocity gradient) are confirmed, and the figure 
also illustrates the significant effect which even very small 
external velocity gradients can have on the velocity on the 
axis. Negative velocity gradients lower the critical swirl value 
and cause type lb flow to break down sooner. Note that a posi- 
tive velocity gradient can prevent breakdown and convert a flow 
from type lb to la, even though the initial profile would 
indicate breakdown in a zero velocity gradient. 

Figure 5.2 presents results of the application of different 
circulation gradients to a type la vortex. The effect of a 
moderate positive circulation gradient is qualitatively the same 
as for a positive velocity gradient. However, the effect 
reverses in type lb flow (Fig. 5.3). Here a positive external 
circulation gradient results in earlier failure, while a nega- 
tive gradient may prevent failure altogether. Note that this 
reversal of the circulation gradient effect can be used to 
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0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 


Fig. 5.2 Development of velocity on the axis as a function 

of external circulation gradient. Initially uniform 
axial flow of type la with Sj[ = 1.346 and 
U e = 1.0 = const. 



X 

Fig. 5.3 Development of velocity on the axis for type lb flow 
under different external axial velocity and circula- 
tion gradients. Initially uniform axial flow. 

U ei = 1.0, K e± = 1.27, Si = 1.425. 

— U e = const, K e = const > U e = variable, K e = const 

— Up = const, K e = variable. 



determine the critical value S Q separating wake-type (la) from 
breakdown- type (lb) behavior. 

5 . 2 Type 2 Vortex Flow 

The effects of external axial velocity and circulation 
gradients on a type 2 vortex with initially uniform axial flow 
are illustrated in Fig. 5.4. For zero gradients the velocity on 
the axis would accelerate, eventually resulting in implosive 
failure, as noted earlier. A positive external axial velocity 
gradient accelerates the increase in velocity on the axis and 
the subsequent failure. A positive velocity gradient applied 
farther downstream has a retarding effect, however. A negative 
external axial velocity gradient initially causes a slower 
increase in velocity on the axis, but subsequently the slope 
begins to steepen, and (implosive) breakdown occurs earlier than 
with zero external velocity gradient. 



Fig. 5.4 Development of velocity on the axis for type 2 flow 
under different external axial velocity and circula- 
tion gradients. Initially uniform axial flow. 

U ei = 1.0, K ei = 3.0, Si = 3.36. 

— U e = const, K e = const; — U e = variable, K e = const; 
U e = const, K e = variable. 
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A positive external circulation gradient hinders the 
acceleration of velocity on the axis. It may even keep the 
velocity on the axis approximately constant and suppress the 
development of extreme velocity peaks in the axial velocity pro- 
file, thus stabilizing the vortex for considerable distance. 

However, in such cases, the velocity on the axis shows a 
divergent oscillatory behavior, and eventual (implosive) failure 
must be expected. Negative external circulation gradients lead 
to further acceleration of the velocity on the axis and to 
earlier failure. 

5 . 3 Type 3 Vortex Flow 

A type 3 vortex with initially uniform axial flow and zero 
external gradient shows a velocity retardation on the axis which 
leads to eventual stagnation there and an inner region of 
reversed axial flow surrounded by an annular region where the 
axial velocity exceeds the freestream velocity. Failure 
eventually results. Figure 5.5 presents the results of external 
axial velocity and circulation gradients on this flow. 

A positive external axial velocity gradient causes further 
deceleration of the velocity on the axis and an earlier failure. 

A negative gradient accelerates the velocity on the axis and 
may retard failure for a significant distance. 

A positive external circulation gradient has the same effect 
as a positive velocity gradient, and it usually leads very 
quickly to a drop in the velocity on the axis and to failure. A 
negative external circulation gradient, if strong enough, will 
accelerate the velocity on the axis and lead to implosive 
failure. 

5 . 4 Type 4 Vortex Flow 

A type 4 vortex with initially uniform axial flow and zero external 
axial velocity or circulation gradients shows increasing velocity 
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Fig. 5.5 Development of velocity on the axis for type 3 flow 
under different external axial velocity and circula- 
tion gradients. Initially uniform axial flow. 

U e± = 1.0, K e . = 7.0, Si = 7.85. 

— U e = const, K e = const; U e = variable, K e = const; 

U e = const, Ke = variable. 

on the axis and, eventually, implosive failure of the inner 
core. The effect of applied external axial velocity and 
circulation gradients is shown in Fig. 5.6. 

A positive external axial velocity gradient initially causes 
slower growth, or deceleration, of the velocity on the axis 
before U again rises steeply and failure occurs. A negative 
external axial velocity gradient, while steepening the rise in 
velocity on the axis, may retard failure for a considerable 
distance . 

A positive external circulation gradient can stabilize 
the flow for some distance. Again, however, a divergent 
oscillation appears in the velocity on the axis, leading to 
eventual failure. A negative external circulation gradient 
leads to quick acceleration of the velocity on the axis and to 
earlier failure. 
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Fig. 5.6 Development of velocity on the axis for type 4 flow 
under different external axial velocity and circula- 
tion gradients. Initially uniform axial flow. 

U ejL - 1.0, K ei = 15.0, Si = 16.81. 

— = const, Ke = const } U e = variable, K e = const 

Ue = const, K e = variable. 

5 . 5 Avoidance of Singularities 

Table 5.1 summarizes the results of subsections 1 through 4. 
Note that the effects of gradients superimpose on the result for 
zero gradients, and that the net result may be an initial 
decrease and later increase in velocity on the axis, for example 
(as in some cases 4 with a positive external axial velocity 
gradient) . 

Material presented here has shown that (1) the core flow is 
mainly responsible for eventual failure, and that (2) the core 
flow is easily influenced by external flow conditions. It is 
therefore obviously possible to influence the development of 
all vortex types by application of external axial velocity or 
circulation gradients. A question of practical importance is 
whether such gradients can be used to stabilize vortex flows 
which are initially of unstable type. The singularity analysis 
of Sec. Ill indicates that this should be possible. 
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Tables, j. 


Behavior of Vortex Plows 


Separating swirl parameter 0 S Q S 2 


Type (uniform initial 
axial flow) 

i 

la lb 

1 

2 

3 

4 


Initial development for zero 
external axial velocity and 
circulation gradients: 
Velocity on axis 

1 

1 

decreases 

1 


decreases 

increases 


Positive external axial 
velocity gradient causes 
additional 

1 

acceleration 

I 

acceleration 

(later: 

deceleration) 

decel . 

deceleration 


(Modest) positive external 
circulation gradient causes 
additional 

i 

i 

accel. decel. 

i 

deceleration 

decel k 

deceleration 



Permanent stabilization of flow initially of type lb by 
external axial velocity or circulation gradients can be achieved, 
as evidenced by Fig. 5.3 and supported by experimental evidence. 
As a further example, Fig. 5.7 demonstrates the result of 
application of a very small positive external axial velocity 
gradient on initially breakdown-unstable type lb vortex flow. 

The gradient was applied at X in the immediate vicinity of 
the failure point, transforming the vortex into a breakdown- 
stable la type. 

It has already been said (Sec. I) that application of the 
parabolic system (1) to subcritical vortex flow (S > S^) is open 
to criticism, but results are nevertheless believed to describe 
trends correctly. The following discussion of behavior of flows 
in the subcritical region is based on the parabolic system and 
refers to behavior of corresponding solutions. 

Since the parabolic system exhibits singularities (S^, t 
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Pig. 5.7 Prevention of breakdown of type lb flow by application 
of a small positive external axial velocity gradient 
beginning at X Q . Initially uniform axial flow, 

Si = 1.402. 

S^/ etc.) it will not be possible to go .from one flow type to 
another without failure of the computation. In particular, it 
will not be possible to go from type 2, 3, or 4 flow to break- 
down-stable type la flow without computational failure. 
Stabilization within a given flow type by application of external 
axial velocity and/or circulation gradients is another matter. 
Failure may be avoided for a considerable distance; and 
permanent stabilization appears possible by application of 
specific external axial velocity and circulation gradients 
(Sec. III). In applications, stable vortex flows are usually 
required only over finite distances , and application of favor- 
able external axial velocity or circulation gradients can result 
in retarding failure for a sufficient distance, as shown in 
previous subsections . 

Stable multi-layer vortex flows of high swirl S > S^ do 
exist in nature (certain tornados and water spouts) and can be 
generated in the laboratory (Rosenzweig, Ross, Lewellen 1962). 
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Their correct and full treatment , beyond establishing general 
trends , appears to require more than a solution of the parabolic 
system (2.1). 
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VI. CONCLUSIONS 


The computation of viscous quasi-cylindrical vortex flows 
at low swirl values does not present any particular dif ficulities 
and can be handled by methods which are only slightly more 
complicated than comparable methods of boundary layer computa- 
tion. Computational difficulties arise at high swirl values 
due to the appearance of swirl- and profile-dependent singulari- 
ties which indicate a failure of the quasi-cylindrical 
approximation and appear to correspond to physical axisymmetric 
• vortex breakdown ' . 

On the basis of recent optical velocity measurements in 
vortex breakdown flows, a distinction has been made in this 
report in particular between 'vortex breakdown', 'vortex jump', 
and the 'vortex bubbles' assocated with vortex breakdown. While 
the vortex breakdown and vortex bubble flow must be computed by 
using the (elliptic) Navier-Stokes equation or the correct 
(inviscid) subset thereof, the flow approaching the breakdown 
can be computed by the (parabolic) quasi-cylindrical vortex 
equations used in this report. This set has been shown to have 
an infinite number of discrete singularities. A 'vortex jump' 
in the sense of Benjamin's (1962) hydraulic jump analogy may 
occur downstream of the vortex breakdown and vortex bubble. A 
solution in the jump region must include a proper matching 
condition between two conjugate flows. 

The present report has presented 

(1) an accurate and efficient method for the computation 
of viscous incompressible quasi-cylindrical axisym- 
metric steady vortex flow, 

(2) a method for obtaining the singular swirl values 
for a given vortex velocity profile combination 
(axial velocity profile and swirl velocity profile, 
or circulation profile) , 
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(3) results of vortex computations for a wide range of 

swirl parameters, initial profiles, and external 

axial velocity and circulation gradients . 

The numerical method of integration of the vortex equations 
is an application of the method of weighted residuals. The 
present method uses exponentials in both the approximating 
expressions for axial velocity and circulation and as weighting 
functions. It is formulated for arbitrary N and can produce 
quick qualitative studies for N = 1 or 2 , or more accurate 
results for N > 2. The results presented were for N = 3. 
Computation of a multitude of vortex flows has been efficient 
and virtually trouble-free. In exploratory investigations such 
as this one and in more detailed analyses of flows, the use of 
N-parameter integral methods appears to have distinct advantages 
and much potential. The singularities of a given vortex velo- 
city profile combination were found by introducing its describing 
parameters into the set of equations for the parameter gradients 
and determining the swirl parameters for which the system 
becomes singular (and the gradients infinite) . 

Major results of the investigation were as follows: 

(1) The equations of viscous quasi-cylindrical incompressible 
vortex flow have an infinite number of discrete singularities. 
These singularities are best categorized in terms of the swirl 
parameter of an equivalent cylinder of fluid in rigid rotation. 
For given vortex velocity profiles, the singular swirl parameter 
values can be computed to determine flow behavior, or results 
such as the present can be used as guidelines . 

(2) The singular swirl parameters S^, , . . . appear to 

correspond to the critical values j^ n /2 of the swirl parameter 
in the solution of the equation of inviscid flow with initial 
rigid rotation. (The regions of validity of the viscous and 
inviscid sets overlap at and near the axis, but the viscous 
quasi-cylindrical set becomes invalid where stream surfaces 


56 



expand or contract rapidly.) The exact values of S^, S^, . . . 
are profile-dependent. 

(3) At each singularity, all derivatives become infinite in 
magnitude and reverse their signs. The result is contrasting 
behavior of vortex flows on both sides of the singularity. In 
particular, as the first singularity is crossed from below 
(increasing swirl parameter) , the tendency to. deceleration of 
the velocity on the axis (and 'expansion of stream surfaces) 
changes to one of acceleration (and contraction) . 

(4) For swirl parameters less than a dividing value S Q < S^, 
vortex flows show smooth viscous decay of all velocity profiles 
to the surrounding freestream velocity (type la) . For swirl 
parameters greater than S Q , increasingly rapid deceleration on 
the axis develops, the swirl paramter S(X) increases further, 
and the computation fails when the singularity S (X) = is 
reached (type lb) . S Q is profile-dependent and appears to 
correspond to the dividing value S Q = /2 for inviscid flow in 
initially rigid rotation. 

(5) External gradients of axial velocity or circulation can 

affect the swirl parameter S (X) and corresponding flow behavior 
by either speeding the approach to a singularity or by avoiding 
it altogether. In particular, breakdown-prone flows of type lb 
with S (X) > S Q can be transformed to type la, which will not 

break down. At the dividing point S Q , the effect of an external 
circulation gradient reverses. 

(6) Each new singularity first enters at the axis and spreads 

outward as S is increased. This leads to alternately decelera- 
ting and accelerating layers in the axial flow profiles. Thus 
for < S < (type 4 flow) , four layers exist: accelerating 

core-flow, adjacent decelerating layer, an accelerating layer 
next, and a decelerating outer layer. It does not appear that 
breakdown-stable vortex configurations are possible under these 
conditions for S > S Q , unless external axial velocity and 
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circulation gradients are applied. 

(7) The singularities are coupled with explosive ot implosive 

expansion, respective contraction of the stream surfaces, and 

the parabolic viscous subset of the Navier-Stokes equations used 

in the present work becomes invalid and should be replaced by 

the inviscid equations of rotating flow or by the full Navier- 

Stokes equations very close to and at the singularities (Bossel 

1969) . These sets permit continuous solutions for all S. The 

results presented here are valid almost to the point of 

computational failure in most cases since stream surface angles 

remain of the order of a few degrees for typical core Reynolds 

4 

numbers of order 10 . 

(8) No qualitative difference exists in the behavior of flows 
of same type but having different velocity profiles (i.e. 
uniform initial axial flow, leading edge vortex, trailing vortex) . 
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APPENDIX 


Coefficients in the ordinary differential equations (2.5) 
N 

= X a * S l,k + U e S 2,k + U ax S 3,k 


n ,k 


n ,k 


B 


n,k 


1=1 

N 


B n,k * X b * Q l,l /2 + K e Q l,2 /2 


1=1 

N 


°k “2 a A S 7 ,k + U e S 8 ,k + U ax S 9 ,] 


1=1 

N 


" X a ^ S 4,k + U e S 5,k + U ax S 6 ,1 


1=1 

N 


= 2 b A,4 /2 + K e Q l, 5 / 2 


i=l 

N 


- X a n T l ,k + a e T 2,k + U ax T 3,k 


m=l 

N 


■ X Vi, 


k + K e S 2,k 


1=1 

N 


“ Y* a A,l + U e Q l,2 + U ax Q l,3 


1=1 

N 


Sc b j^7. 


k + K e S 8,k 


£=1 

N 

X 

1=1 


b A S 4,k + K e S 5,k 
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N 


- X a *°l,4 + “e°l,5 + “ax 0 !,. 


£=1 

N 


“ X Vi, 


i + k t_ ; 
k- e 2,k 


n=l 


These coefficients involve present values of the parameters 
a n' ^n' an< ^ °f u _ / U o , and K^. They must therefore be calcula- 


ax 


ted anew at each step. However, the numbers Q, S, T, S, and T 
are constants which are determined only once at the beginning 
of the computation. It is convenient to define them in terms 
of other constants Z, P, and R. The following numbers are all 
for given £ and n. The index k runs from 1 to N + 1. 

Let 


z l,l 

= 

a k : 

J- 


Z 2,l 

= 

(a k 

+ 

c,)- 1 

Z 3 , 1 

= 

(a k 

+ 

2a) 1 

Z 4 , 1 

= 

(a k 

+ 

na) 1 

Z 5,l 

= 

[ °k 

+ 

(n + Do] -1 

Z 6,l 

= 

“k 

+ 

(n + 2) a]" 1 

Z 7,l 

= 

[ *k 

+ 

(n + £) a] 1 

Z 8,l 

= 

[ *k 

+ 

(n + £ + l)a] -1 

Z 9,l 

= 

<‘k 

+ 

(n + £ + 2) a] 1 

Z 10,l 

= 

[a k 

+ 

£a] _1 

Z ll,l 

= 


+ 

(£ + l)a] _1 
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Z 12,l = 

[a k + U + 2) a] 1 


Z i,2 = 

Z. n Z. . 
i/l i/l 



Z i , 3 = 

2Z. „Z. n 
1,2 i,l 

► i = 1, 2 , . . . , 


Z . . 
1/4 

3Z . -Z. i 
1,3 1,1 


Then 


1 



P m, 1 " 

z„ - z r 
4 ,m 5,m 



P m, 2 

Z l,m Z 2,m 

► m = 1, 2, 3, 4 


P m, 3 

Z 2,m 


For the Q ' s 

and R's the 

index m = 1 , 2 

/ 3: 

, 1 

= Z 7 ,m 

2Z S + Z Q 

8 f in S? f m 


®m, 2 

= Z 4,m " 

^ Z 5,m + Z 6,m 


^m, 3 

Z 5 ,m 

Z 6 ,m 


4 

Z 10 ,m 

Z-, 

11 ,m 


^m,5 

= Z 1 ,m 

2Z_ + Z- 

2 ,m 3 ,m 


Q m,6 

~ Z 2 ,m 

Z-, 

3 ,m 


®m,7 

= Z 11 ,m “ 

Z 12,m 


00 

e 

a 

~ Z 2 ,m 

z, 

3 ,m 


®m,9 

= Z 3,m 



R m, 1 

- ^ Z 10,m 

' Z 8,m - Z ll,m + Z 9,m ) '' l < n + 


+ ( - z io 

_ + z_ + z nl 

f m 7 fin 11 j 

- Z_ ) / (na) 

rin 8 f m 


V • 

: 
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Then 


8* 

lO 

= 

< Z l,m - 

Z. m ~ Z 9 m + J / t < n + D “3 

5 /in. 2 ,m d ,m 



+ (-z x< 

m + Z 4,m + Z 2,m " Z 5,m )/(na) ' 

R m,3 

S 

(Z 2,m - 

z 6 /m )/[ (n + l)a] + (-Z 2fm + Z 5>m )/(na) 

R m,4 

= 

(Z 10,m 

2Z ll,m + Z 12,m )/a ” Z 10,m+1 + Z ll,m+1 

R m, 5 

= 

(Z, m " 
1 ,m 

2Z 2 ,m + Z 3,m )/a " Z l,m+1 " Z 2,m+1 

R m,6 

= 

(Z 2,m " 

Z 3,m )a ” Z 2,m+1 

R m, 7 

= 

( ” Z 10,rn 

+ 2Z, , - Z, 0 _)/a 

H f in o. 2 

R m, 8 

s= 

^” Z l,m 

+ 2Z 2 ,m ” Z 3,m )/o ‘ 

R m,9 

= 

(Z 3,m " 

Z 2,m )/a 

S i,k 

= 

2a k Q 3,i 

4Q 2,i + a k R 3 ,i “ 4a k R 2,i + 2R l,i 

§ i/k 

= 

Q m . + 

l/l 

^ 1=1,2,. . .9 


T i,k “ - 2o k 3p 4,i + 14o k 2p 3,i ' 20o k P 2,i + 4P l,i 


\,k ‘ - 2o k P 2,i + 4a k P l,i 


i = 1, 2, 3 
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